## *************************************
##          B2 More Controls       
##        Hao Zhang & Ye Zhang
##             2025/7/24
## *************************************

# load packages
library(lfe)
library(stargazer)
library(tidyverse)

# load data
rm(list = ls())
load("../analysis data/firm analysis.RData")

# run inverse hyperbolic sine
ihs <- function(x) {
  y <- log(x + sqrt(x ^ 2 + 1))
  return(y)
}

# unemployment insurance
m1 <- felm(insurance_dummy ~ gdp_pressure + ihs(unemp_rate) + ihs(employment) + 
             ihs(wage) + ihs(output1) + ihs(profit) +  
             log(FDI + 1) + gdpgrowth + log(gdp/population) +
             ihs(deficit/gdp) +  ihs(capacity) + tax  | 
             frdm + year + educ + ethnicity  + gender + 
             begin + connection + seq | 0 | citycode, data = firm)

firm1 <- firm %>% filter(rate1 > 0)
m2 <- felm(log(rate1) ~ gdp_pressure + ihs(unemp_rate) + ihs(employment) +
             ihs(wage) + ihs(output1) + ihs(profit) +  
             log(FDI + 1) + gdpgrowth + log(gdp/population) +
             ihs(deficit/gdp) +  ihs(capacity) + tax  | 
             frdm + year + educ + ethnicity  + gender + 
             begin + connection + seq | 0 | citycode, data = firm1)

# medical insurance
m3 <- felm(medical_dummy ~ gdp_pressure + log(hospital/population) + ihs(employment) +
             ihs(wage) + ihs(output1) + ihs(profit) +  
             log(FDI + 1) + gdpgrowth + log(gdp/population) +
             ihs(deficit/gdp) +  ihs(capacity) + tax | 
             frdm + year + educ + ethnicity  + gender + 
             begin + connection + seq | 0 | citycode, data = firm)

firm2 <- firm %>% filter(rate2 > 0)
m4 <- felm(log(rate2) ~ gdp_pressure + log(hospital/population) + ihs(employment) +
             ihs(wage) + ihs(output1) + ihs(profit) +  
             log(FDI + 1) + gdpgrowth + log(gdp/population) +
             ihs(deficit/gdp) +  ihs(capacity) + tax | 
             frdm + year + educ + ethnicity  + gender + 
             begin + connection + seq | 0 | citycode, data = firm2)

# output regression table
stargazer(m1, m2, m3, m4, title = "Interjurisdictional Competition and Social Insurances (More Controls)", 
          dep.var.labels = c("Unemployment", "Unemployment", "Health and Pension", "Health and Pension"),
          omit.stat = c("rsq", "adj.rsq", "ser"), 
          font.size = "small", 
          add.lines = list(c("Mayor Characteristics", "Y", "Y", "Y", "Y"),
                           c("Firm Characteristics", "Y", "Y", "Y", "Y"), 
                           c("Firm Fixed Effects", "Y", "Y", "Y", "Y"), 
                           c("Year Fixed Effects", "Y", "Y", "Y", "Y")))
